home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Programming / Source / HippoDraw / hippo / hippoplot.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-04-28  |  56.9 KB  |  2,627 lines
  1. /*
  2.  * hippoplot.c - Routines for displaying hippo ntuples.
  3.  *
  4.  * Copyright (C)  1991  The Board of Trustees of The Leland Stanford
  5.  * Junior University.  All Rights Reserved.
  6.  *
  7.  * $Id: hippoplot.c,v 3.20 1992/04/24 01:12:19 rensing Rel $
  8.  *
  9.  * by jonas karlsson, at SLAC, August 1990
  10.  *  split up by Paul Rensing, Feb 28,1991
  11.  *  modified by M. Gravina, March 28, 1991
  12.  */
  13.  
  14. #include <stdio.h>
  15. #include <stdlib.h>
  16. #include <math.h>
  17. #include <string.h>
  18. #include <float.h>
  19. #include <ctype.h>
  20. #include "hippo.h"
  21. #include "hippoutil.h"
  22.  
  23. GLOB_QUAL const char hippoplot_c_rcsid[] = 
  24.      "$Id: hippoplot.c,v 3.20 1992/04/24 01:12:19 rensing Rel $";
  25.  
  26.  
  27. #ifdef _NEXT_PLOT_
  28. #include "hippoplotNeXT.h"
  29. #ifndef DEF_PLOT_DRVR
  30. #define DEF_PLOT_DRVR NEXT
  31. #endif
  32. #endif
  33.  
  34. #ifdef _UNIXPLOT_PLOT_
  35. #include "hippoplotUP.h"
  36. #ifndef DEF_PLOT_DRVR
  37. #define DEF_PLOT_DRVR UNIXPLOT
  38. #endif
  39. #endif
  40.  
  41. #ifdef _XIV_PLOT_
  42. #include "hippoplotXIV.h"
  43. #ifndef DEF_PLOT_DRVR
  44. #define DEF_PLOT_DRVR XIVPLOT
  45. #endif
  46. #endif
  47.  
  48. #ifdef _X11_PLOT_
  49. #include "hippoplotX11.h"
  50. #ifndef DEF_PLOT_DRVR
  51. #define DEF_PLOT_DRVR X11PLOT
  52. #endif
  53. #endif
  54.  
  55. #ifdef THINK_C
  56. #include "hippoplotMAC.h"
  57. #ifndef DEF_PLOT_DRVR
  58. #define DEF_PLOT_DRVR MAC
  59. #endif
  60. #endif
  61.  
  62. #include "hippoplotPS.h"
  63. #ifndef DEF_PLOT_DRVR
  64. #define DEF_PLOT_DRVR PSPLOT
  65. #endif
  66.  
  67. #define NUM_FUZZ FLT_EPSILON*4
  68.  
  69. /* find offset to matrix element */
  70. #define INDEX(row, column, totcols) ((row) * (totcols) + (column))
  71. #define BININDEX(x, y, z, xs, ys) ((xs) * (ys) * (z) + (x) * (ys) + y)
  72.  
  73. #define MIN(x, y)  (((x) > (y)) ? (y) : (x))
  74. #define MAX(x, y)  (((x) > (y)) ? (x) : (y))
  75.  
  76. #define XPADDING(disp) ((disp)->drawRect.size.width*0.01)
  77. #define YPADDING(disp) ((disp)->drawRect.size.height*0.01)
  78.  
  79.  
  80. #ifdef VMS
  81. #define Concat3(x,y,z) x/**/y/**/z
  82. #else
  83. #define Concat3(x,y,z) x##y##z
  84. #endif
  85.  
  86. #ifdef _NEXT_PLOT_
  87. #define NeXTFunc(name,parm)           \
  88. case NEXT:                            \
  89.      rc = Concat3(name,_NeXT,parm);   \
  90.      break;
  91.  
  92. #else
  93. #define NeXTFunc(name,parm)  
  94. #endif
  95.  
  96. #ifdef _UNIXPLOT_PLOT_
  97. #define UPFunc(name,parm)           \
  98. case UNIXPLOT:                      \
  99.      rc = Concat3(name,_UP,parm);   \
  100.      break;
  101.  
  102. #else
  103. #define UPFunc(name,parm)  
  104. #endif
  105.  
  106. #ifdef _XIV_PLOT_
  107. #define XIVFunc(name,parm)         \
  108. case XIVPLOT:                      \
  109.      Concat3(name,_XIV,parm);      \
  110.      rc = 0;                       \
  111.      break;
  112.  
  113. #else
  114. #define XIVFunc(name,parm)  
  115. #endif
  116.  
  117. #ifdef _X11_PLOT_
  118. #define X11Func(name,parm)           \
  119. case X11PLOT:                        \
  120.      Concat3(name,_X11,parm);        \
  121.      rc = 0;                         \
  122.      break;
  123.  
  124. #else
  125. #define X11Func(name,parm)  
  126. #endif
  127.  
  128. #ifdef THINK_C
  129. #define MACFunc(name,parm)           \
  130. case MAC:                            \
  131.      Concat3(name,_MAC,parm);        \
  132.      rc = 0;                         \
  133.      break;
  134.  
  135. #else
  136. #define MACFunc(name,parm)  
  137. #endif
  138.  
  139. /*
  140.  * always make PS driver available
  141.  */
  142. #define PSFunc(name,parm)           \
  143. case PSPLOT:                        \
  144. case EPSPLOT:                       \
  145.      Concat3(name,_PS,parm);        \
  146.      rc = 0;                        \
  147.      break;
  148.  
  149. #define PlotSwitch(drvr,name, parm) \
  150. switch (drvr)                       \
  151. {                                   \
  152.      NeXTFunc(name,parm)            \
  153.      UPFunc(name,parm)              \
  154.      XIVFunc(name,parm)             \
  155.      X11Func(name,parm)             \
  156.      MACFunc(name,parm)             \
  157.      PSFunc(name,parm)              \
  158.      default:                       \
  159.       rc = -1;                  \
  160. }
  161.  
  162. struct maxs 
  163. {
  164.      float xl,xh,yl,yh;
  165. };
  166.  
  167.  
  168. /* 
  169.  * private functions
  170.  */
  171. static float calcTicks(float size, float *magnitude);
  172.  
  173. static int drawData( display disp );
  174. static int drawData1D( display disp );
  175. static int drawAxis( display disp );
  176. static int drawLabels( display disp );
  177.  
  178. static int drawLego1D(display disp, float *lego, int *over, int npts);
  179. static int drawPoint1D(display disp, float *xy, int *over, int npts);
  180. static int drawError1D(display disp, float *xy, float *err, 
  181.                int *over, int npts, int err_type);
  182. static int drawLine1D(display disp, float *xy, int *over, 
  183.               int npts, int xmin_pt, linestyle_t ls);
  184. static void box_cross(float *p1, float *p2, display disp);
  185. int intcmp(const void *i1, const void *i2);
  186.  
  187. static int drawColor2D(display disp);
  188. static int drawScatter2D(display disp);
  189. static int drawLego2D(display disp);
  190.  
  191. static int initPlot(display disp);
  192. static int endPlot(display disp);
  193. static int drawAxisBox( display disp );
  194. static int drawTicks(display disp, binding_t axis);
  195. static int drawXTicks(float *ticks, int nticks, float tickLength, char tickLoc );
  196. static int drawYTicks(float *ticks, int nticks, float tickLength, char tickLoc );
  197. static int drawMag( float x, float y, int mag, float fontSize );
  198.  
  199. static int drawText(char *s, float x, float y, float fontsize, float angle,
  200.             char xp, char yp );
  201.  
  202. static int getXYData( display disp, float **xy, float **xerr, 
  203.     float **yerr, int **over, struct maxs *m, int *xmin_pt);
  204. static int getHistoData( display disp, float **xy, float **xylego, 
  205.     float **xerr, float **yerr, int **over,struct maxs *m);
  206. static int getScatData( display disp, float **xy, struct maxs *m );
  207.  
  208. static int plotFunc(display disp, func_id fun );
  209.  
  210. /*
  211.  * global variables.
  212.  */
  213. static int init_done;
  214. static plotdrvr_t plot_drvr = DEF_PLOT_DRVR;
  215.  
  216. #ifdef _XIV_PLOT_
  217. /* For passing to the X Interviews driver... */
  218. static void* canvas;
  219. static void* painter;
  220. #endif
  221.  
  222. #ifdef _X11_PLOT_
  223. /* for passing to X11 driver ... */
  224. static void* dpy;
  225. static void* screen;
  226. static void* drawable;
  227. static void* gc;
  228. #endif
  229.  
  230.  
  231. int h_setPlotDrvr( plotdrvr_t drvr, ...)
  232. {
  233.      int rc;
  234.      va_list argPtr;
  235.      void *parm;
  236.      
  237.      switch (drvr)
  238.      {
  239. #ifdef _NEXT_PLOT_
  240.      case NEXT:
  241. #endif
  242. #ifdef _XIV_PLOT_
  243.      case XIVPLOT:
  244. #endif      
  245. #ifdef _X11_PLOT_
  246.      case X11PLOT:
  247. #endif      
  248. #ifdef THINK_C
  249.      case MAC:
  250. #endif
  251.      case LPR:
  252.       rc = 0;
  253.       break;
  254.             
  255.      case PSPLOT:
  256.      case EPSPLOT:
  257.       va_start( argPtr, drvr );
  258.       parm = va_arg(argPtr, FILE *);
  259.       rc = initDrvr_PS( parm );
  260.       va_end(argPtr);
  261.       break;
  262.  
  263. #ifdef _UNIXPLOT_PLOT_
  264.      case UNIXPLOT:
  265.       va_start( argPtr, drvr );
  266.       parm = va_arg(argPtr, FILE *); /* this shit does not work yet */
  267.       rc = initDrvr_UP( parm );
  268.       va_end(argPtr);
  269.       break;
  270. #endif
  271.             
  272.      default:
  273.       h_error("h_setPlotDrvr - invalid plot driver. Driver may not be compiled in this version.\n");
  274.       rc = -1;
  275.      }
  276.      
  277.      if (rc == 0) 
  278.       plot_drvr = drvr;
  279.      else
  280.       h_error("h_setPlotDrvr - error in driver init. Driver not changed.\n");
  281.      
  282.      return rc;
  283. }
  284.  
  285. int h_endPlotDrvr( void )
  286. {
  287.      int rc;
  288.      
  289.      switch (plot_drvr)
  290.      {
  291.      case PSPLOT:
  292.      case EPSPLOT:
  293.       rc = endDrvr_PS( );
  294.       break;
  295.  
  296.      default:
  297.       rc = 0;
  298.      }
  299.      
  300.      return rc;
  301. }
  302.  
  303. int h_endPage( void )
  304. {
  305.      int rc;
  306.      
  307.      switch (plot_drvr)
  308.      {
  309.      case PSPLOT:
  310.      case EPSPLOT:
  311.       rc = endPage_PS( );
  312.       break;
  313.  
  314.      default:
  315.       rc = 0;
  316.      }
  317.      
  318.      return rc;
  319. }
  320.  
  321. int h_plot(display disp,...)
  322. { 
  323.      int rc = 0;
  324.      func_id fun;
  325. #if defined(_XIV_PLOT_) || defined(_X11_PLOT_)
  326.      va_list argPtr;
  327. #endif
  328.      
  329.      /* 
  330.       * be nice. If driver is lineprinter, call h_print.
  331.       */
  332.      if (plot_drvr == LPR)
  333.      {
  334.       h_print( disp );
  335.       return 0;
  336.      }
  337.      
  338. #ifdef _XIV_PLOT_
  339.      if (plot_drvr == XIVPLOT)
  340.      {
  341.       va_start( argPtr, disp );
  342.       painter = va_arg(argPtr, void *);
  343.       canvas  = va_arg(argPtr, void *);
  344.       va_end(argPtr);
  345.      }
  346. #endif
  347.  
  348. #ifdef _X11_PLOT_
  349.      if (plot_drvr == X11PLOT)
  350.      {
  351.       va_start( argPtr, disp );
  352.       dpy       = va_arg(argPtr, void *);
  353.       screen    = va_arg(argPtr, void *);
  354.       drawable  = va_arg(argPtr, void *);
  355.       gc        = va_arg(argPtr, void *);
  356.       va_end(argPtr);
  357.      }
  358. #endif
  359.  
  360.      if (h_bin(disp) != 0) return -1;
  361.      
  362.      init_done = 0;
  363.           
  364.      drawData(disp);
  365.      
  366.      for (fun = disp->plot_func; fun != NULL && rc == 0; fun = fun->next )
  367.      {
  368.       rc = plotFunc( disp, fun );
  369.      }
  370.  
  371.      if (disp->flags.drawAxes) drawAxis(disp);    
  372.      
  373.      if (disp->flags.drawTitles) drawLabels( disp );
  374.      
  375.      endPlot( disp );
  376.      
  377.      return 0;
  378. }
  379.  
  380.  
  381. #define LABEL_SPACE_X    4.0
  382. #define LABEL_SPACE_Y    4.0
  383. #define MIN_TICKS 4
  384. #define MAX_TICKS (MIN_TICKS*2+1)
  385.  
  386. static float calcTicks(float size, float *magnitude)
  387. {
  388.      static float    goodTicks[] = {10.0, 5.0, 4.0, 2.0, 1.0};
  389.      float        tickSize;
  390.      int        tickIndex;
  391.  
  392.      if (size <= 0.0) 
  393.      {
  394.       h_error("calcTicks: size is <= 0");
  395.       size = fabs(size);
  396.       if (size == 0.0) size = 1.0;
  397.      }
  398.      
  399.      *magnitude = floor(log10(size));
  400.      if (size/pow(10.0,*magnitude) < MIN_TICKS) (*magnitude)--;
  401.  
  402.      /*
  403.       *  now fit the max number of ticks into this range
  404.       */
  405.      for (tickIndex = 0;
  406.       size / (tickSize=goodTicks[tickIndex]*pow(10.0,*magnitude) ) 
  407.       < MIN_TICKS;
  408.       tickIndex++);
  409.  
  410.      if (tickIndex == 0) 
  411.       (*magnitude)++;
  412.      
  413.      return tickSize;
  414. }
  415.  
  416. static int drawAxis(display disp)
  417. {
  418.      if (disp == NULL) return -1;
  419.      
  420.      if (disp->graphtype == LEGOPLOT && disp->dim == 2) return 0;
  421.      
  422.      if (!init_done)
  423.      {
  424.       initPlot(disp);
  425.       init_done = 1;
  426.      }
  427.      
  428.      drawAxisBox( disp );
  429.      if (drawTicks(disp,YAXIS))
  430.       h_error("drawAxis: error drawing y ticks");
  431.      if (drawTicks(disp,XAXIS)) 
  432.       h_error("drawAxis: error drawing x ticks");
  433.  
  434.      return 0;
  435. }
  436.  
  437.  
  438. static int drawData(display disp)
  439. {
  440.      float    fontHeight;
  441.  
  442.      if (disp->ntuple == NULL) return 0;
  443.  
  444.      if (disp->ntuple->ndata == 0)
  445.      { 
  446.       fontHeight= MIN(disp->drawRect.size.width*0.05,24.0);
  447.       drawText("No Data in Tuple",
  448.            disp->marginRect.origin.x+disp->marginRect.size.width/2.0,
  449.            disp->marginRect.origin.y+disp->marginRect.size.height/2.0,
  450.            fontHeight, 0, 'c', 'c');
  451.       return 0;
  452.      }
  453.      
  454.      switch (disp->graphtype) 
  455.      {
  456.      case HISTOGRAM:
  457.      case XYPLOT:
  458.      case STRIPCHART:
  459.       drawData1D( disp );
  460.       break;
  461.  
  462.      case COLORPLOT:        /* 2D COLOR */
  463.       drawColor2D(disp);
  464.       break;
  465.       
  466.      case SCATTERPLOT:        /* 2D SCATTER */
  467.       drawScatter2D(disp);
  468.       break;
  469.       
  470.      case LEGOPLOT:        /* 2D LEGO */
  471.       drawLego2D(disp);
  472.       break;
  473.      }
  474.      
  475.      
  476.      return 0;
  477. }
  478.  
  479. static int drawData1D(display disp)
  480. {
  481.      float *xy = NULL, *lego = NULL, *xerr = NULL, *yerr = NULL;
  482.      int *over = NULL, xmin_pt = 0, npts;
  483.      struct maxs limits;
  484.      
  485.      /*
  486.       * first get the appropriate data, depending on the graphtype.
  487.       */
  488.      switch (disp->graphtype) 
  489.      {
  490.      case HISTOGRAM:        /* Histogra */
  491.       npts = getHistoData(disp,&xy,&lego,&xerr,&yerr,&over,&limits);
  492.       if (npts == 0) return 0;
  493.       if (npts < 0) 
  494.       {
  495.            h_error("h_plot - error getting bin data");
  496.            return -1;
  497.       }
  498.  
  499.       if (disp->yAxis.flags.autoScale)
  500.       {
  501.            if (limits.yh < limits.yl)
  502.            {
  503.             /* must be no points in range. */
  504.             return 0;
  505.            }             
  506.            else if (limits.yh == limits.yl)
  507.            {
  508.             /* must be only on point in range. */
  509.             disp->yAxis.low = limits.yh - 0.5;
  510.             disp->yAxis.high = limits.yh + 0.5;
  511.            }
  512.            else
  513.            {
  514.             if (limits.yh > 0.0)
  515.              disp->yAxis.high = limits.yh;
  516.             else
  517.             {
  518.              if (disp->yAxis.flags.log)
  519.              {
  520.                   h_error("All data < 0 on log axis");
  521.                   return -1;
  522.              }
  523.              disp->yAxis.high = 0.0;
  524.             }
  525.             
  526.             if (limits.yl < 0.0 || disp->yAxis.flags.log)
  527.              disp->yAxis.low = limits.yl;
  528.             else
  529.              disp->yAxis.low = 0.0;
  530.            }
  531.       }
  532.       else if (disp->yAxis.flags.log && disp->yAxis.low <= 0.0)
  533.       {
  534.            if (limits.yl > 0.0)
  535.             disp->yAxis.low = limits.yl;
  536.            else
  537.             disp->yAxis.low = pow(10.0, FLT_MIN);
  538.            if (disp->yAxis.high <= disp->yAxis.low)
  539.             disp->yAxis.high = disp->yAxis.low*(1.0+NUM_FUZZ);
  540.       }
  541.  
  542.       break;
  543.       
  544.      case XYPLOT:
  545.      case STRIPCHART:
  546.       npts = getXYData(disp,&xy,&xerr,&yerr,&over,&limits,&xmin_pt);
  547.       if (npts == 0) return 0;
  548.       if (npts < 0) 
  549.       {
  550.            h_error("h_plot - error getting data for x-y plot");
  551.            return -1;
  552.       }
  553.  
  554.       if (disp->graphtype == XYPLOT) xmin_pt = 0;
  555.       
  556.       /* 
  557.        * autoscale axes 
  558.        */
  559.       if (disp->xAxis.flags.autoScale)
  560.       {
  561.            if (limits.xh > limits.xl)
  562.            {
  563.             disp->xAxis.low = limits.xl;
  564.             disp->xAxis.high = limits.xh;
  565.            }
  566.            else if (limits.xh == limits.xl)
  567.            {
  568.             /* must be only on point in range. */
  569.             disp->xAxis.low = limits.xl - 0.5;
  570.             disp->xAxis.high = limits.xh + 0.5;
  571.            }
  572.            else
  573.            {
  574.             /* must be no points in range. */
  575.             return 0;
  576.            }             
  577.            h_adjustAxis(&(disp->xAxis.low), &(disp->xAxis.high),
  578.                 0, disp->xAxis.flags.log);
  579.       }
  580.       else if (disp->xAxis.flags.log && disp->xAxis.low <= 0.0)
  581.       {
  582.            if (limits.xl > 0.0)
  583.             disp->xAxis.low = limits.xl;
  584.            else
  585.             disp->xAxis.low = pow(10.0, FLT_MIN);
  586.            if (disp->xAxis.high <= disp->xAxis.low)
  587.             disp->xAxis.high = disp->xAxis.low*(1.0+NUM_FUZZ);
  588.       }
  589.  
  590.       if (disp->yAxis.flags.autoScale)
  591.       {
  592.            if (limits.yh > limits.yl)
  593.            {
  594.             disp->yAxis.low = limits.yl;
  595.             disp->yAxis.high = limits.yh;
  596.            }
  597.            else if (limits.yh == limits.yl)
  598.            {
  599.             /* must be only on point in range. */
  600.             disp->yAxis.low = limits.yl - 0.5;
  601.             disp->yAxis.high = limits.yh + 0.5;
  602.            }
  603.            else
  604.            {
  605.             /* must be no points in range. */
  606.             return 0;
  607.            }             
  608.       }
  609.       else if (disp->yAxis.flags.log && disp->yAxis.low <= 0.0)
  610.       {
  611.            if (limits.yl > 0.0)
  612.             disp->yAxis.low = limits.yl;
  613.            else
  614.             disp->yAxis.low = pow(10.0, FLT_MIN);
  615.            if (disp->yAxis.high <= disp->yAxis.low)
  616.             disp->yAxis.high = disp->yAxis.low*(1.0+NUM_FUZZ);
  617.       }
  618.       
  619.       break;
  620.       
  621.      default:
  622.       h_error("hippo error: invalid disp->graphtype\n");
  623.       return -1;
  624.       break;
  625.      }
  626.  
  627.      if (disp->yAxis.flags.autoScale)
  628.       h_adjustAxis(&(disp->yAxis.low), &(disp->yAxis.high),
  629.                0, disp->yAxis.flags.log );
  630.      
  631.      /*
  632.       * init the plot device
  633.       */
  634.      if (!init_done)
  635.      {
  636.       initPlot( disp );
  637.       init_done = 1;
  638.      }
  639.  
  640.      /*
  641.       * draw the points, lines etc.
  642.       */
  643.  
  644.      /* draw points */
  645.      if ( disp->drawtype & POINT )
  646.       drawPoint1D(disp, xy, over, npts);
  647.      
  648.      /* draw error bars */
  649.      if ( disp->drawtype & ERRBAR )
  650.      {
  651.       if (yerr != NULL) drawError1D(disp,xy,yerr,over,npts,YERROR);
  652.       if (xerr != NULL) drawError1D(disp,xy,xerr,over,npts,XERROR);
  653.      }
  654.      
  655.      
  656.      /* join points */
  657.      if ( disp->drawtype & LINE )
  658.       drawLine1D(disp,xy,over,npts,xmin_pt,disp->lineStyle);
  659.      
  660.      /* draw bars */
  661.      if ( (disp->drawtype & BOX) && disp->graphtype==HISTOGRAM )
  662.       drawLego1D(disp, lego, over, npts);
  663.      
  664.      free(lego);
  665.      free(xy);
  666.      free(xerr);
  667.      free(yerr);
  668.      free(over);
  669.      
  670.      return 0;
  671. }
  672.  
  673.  
  674. static int initPlot(display disp)
  675. {
  676.      int rc;
  677.      rectangle userRect;
  678.  
  679.      if (disp->xAxis.flags.log)
  680.      {
  681.       userRect.origin.x = log10(disp->xAxis.low);
  682.       userRect.size.width = log10(disp->xAxis.high) - 
  683.            log10(disp->xAxis.low);
  684.      }
  685.      else
  686.      {
  687.       userRect.origin.x = disp->xAxis.low;
  688.       userRect.size.width = disp->xAxis.high - disp->xAxis.low;
  689.      }
  690.      if (disp->yAxis.flags.log)
  691.      {
  692.       userRect.origin.y = log10(disp->yAxis.low);
  693.       userRect.size.height = log10(disp->yAxis.high)
  694.            - log10(disp->yAxis.low);
  695.      }
  696.      else
  697.      {
  698.       userRect.origin.y = disp->yAxis.low;
  699.       userRect.size.height = disp->yAxis.high - disp->yAxis.low;
  700.      }
  701.      
  702.      switch (plot_drvr)
  703.      {
  704.       NeXTFunc(initPlot,(&disp->drawRect, &disp->marginRect, &userRect));
  705.       UPFunc(initPlot,(&disp->drawRect, &disp->marginRect, &userRect));
  706.       MACFunc(initPlot,(&disp->drawRect, &disp->marginRect, &userRect));
  707.       PSFunc(initPlot,(&disp->drawRect, &disp->marginRect, &userRect));
  708.       
  709. #ifdef _XIV_PLOT_
  710.      case XIVPLOT:
  711.       initPlot_XIV(painter,canvas);
  712.       setHistoCoords_XIV(&disp->marginRect, &userRect);
  713.       rc = 0;
  714.       break;
  715. #endif
  716.  
  717. #ifdef _X11_PLOT_
  718.      case X11PLOT:
  719.       initPlot_X11(dpy,screen,drawable,gc);
  720.       setHistoCoords_X11(&disp->drawRect, &disp->marginRect, &userRect);
  721.       rc = 0;
  722.       break;
  723. #endif
  724.  
  725.      default:
  726.       rc = -1;
  727.      }
  728.      
  729.      return rc;
  730. }
  731.  
  732.  
  733.  
  734. static int endPlot(display disp)
  735. {
  736.      int rc;
  737.      
  738.      switch (plot_drvr)
  739.      {
  740.       UPFunc(endPlot,());
  741.             
  742.      default:
  743.       rc = 0;
  744.      }
  745.      
  746.      return rc;
  747. }
  748.  
  749.  
  750.  
  751. static int drawLabels( display disp )
  752. {
  753.      float x, y;
  754.      float fontSize;
  755.      char string[80];
  756.      
  757.      /*
  758.       * title
  759.       */
  760.      fontSize = disp->drawRect.size.width*0.05;
  761.      fontSize = MIN(fontSize, 
  762.             (disp->drawRect.size.height - 
  763.              disp->marginRect.size.height - 2*YPADDING(disp) -
  764.              (disp->marginRect.origin.y-disp->drawRect.origin.y)));
  765.      x = disp->marginRect.origin.x + 0.5 * disp->marginRect.size.width;
  766.      y = disp->drawRect.origin.y + disp->drawRect.size.height - YPADDING(disp);
  767.      h_expandLabel(string,disp->title,80,disp);
  768.      if (fontSize > disp->marginRect.size.width/strlen(string)*2)
  769.      {
  770.       fontSize = disp->drawRect.size.width/strlen(string)*2;
  771.       x = disp->drawRect.origin.x + 0.5 * disp->drawRect.size.width;
  772.      }
  773.      
  774.      if (drawText(string, x, y, fontSize,0.0,'c','t') != 0)
  775.      {
  776.       h_error("Error drawing plot title");
  777.       return -1;
  778.      }
  779.  
  780.      /*
  781.       * x axis label
  782.       */
  783.      fontSize = MIN(disp->drawRect.size.width*0.04,14.0);
  784.      fontSize = MIN(fontSize, 
  785.              (disp->marginRect.origin.y-disp->drawRect.origin.y)/2.0);
  786.      fontSize = MIN(fontSize, 
  787.             disp->drawRect.size.width/strlen(disp->xAxis.label)*2);
  788.      y = disp->drawRect.origin.y + YPADDING(disp);
  789.      h_expandLabel(string,disp->xAxis.label,80,disp);
  790.      if (drawText(string, x, y, fontSize,0.0,'c','b') != 0)
  791.      {
  792.       h_error("Error drawing x axis label");
  793.       return -1;
  794.      }
  795.      
  796.      /*
  797.       * y axis label
  798.       */
  799.      fontSize = MIN(disp->drawRect.size.width*0.04,14.0);
  800.      fontSize = MIN(fontSize, 
  801.             disp->drawRect.size.height/strlen(disp->yAxis.label)*2);
  802.      x = disp->drawRect.origin.x + XPADDING(disp);
  803.      y = disp->marginRect.origin.y + disp->marginRect.size.height/2.0;
  804.      h_expandLabel(string,disp->yAxis.label,80,disp);
  805.      if (drawText(string, x, y, fontSize,90.0,'c','t') != 0)
  806.      {
  807.       h_error("Error drawing y axis label");
  808.       return -1;
  809.      }
  810.      
  811.      
  812.      return 0;
  813. }
  814.  
  815.      
  816.  
  817. static int drawText(char *s, float x, float y, float fontsize, float angle,
  818.          char xp, char yp )
  819. {
  820.      int rc;
  821.      
  822.      PlotSwitch(plot_drvr,drawText,(s, x, y, fontsize, angle, xp, yp));
  823.      
  824.      return rc;
  825. }
  826.  
  827.  
  828. static int drawAxisBox( display disp )
  829. {
  830.      int rc = 0;
  831.      float xy[10];
  832.  
  833.      if (disp->xAxis.flags.log)
  834.      {
  835.       xy[0] = xy[2] = xy[8] = log10(disp->xAxis.low);
  836.       xy[4] = xy[6] = log10(disp->xAxis.high);
  837.      }
  838.      else
  839.      {
  840.       xy[0] = xy[2] = xy[8] = disp->xAxis.low;
  841.       xy[4] = xy[6] = disp->xAxis.high;
  842.      }
  843.      if (disp->yAxis.flags.log)
  844.      {
  845.       xy[1] = xy[7] = xy[9] = log10(disp->yAxis.low);
  846.       xy[3] = xy[5] = log10(disp->yAxis.high);
  847.      }
  848.      else
  849.      {
  850.       xy[1] = xy[7] = xy[9] = disp->yAxis.low;
  851.       xy[3] = xy[5] = disp->yAxis.high;
  852.      }
  853.      
  854.      PlotSwitch(plot_drvr,drawLine,(xy,5, SOLID));
  855.  
  856.      return rc;
  857. }
  858.  
  859. #define N__TICKS 20
  860.  
  861. static int drawTicks(display disp, binding_t axis)
  862. {
  863.      int nticks = 0;
  864.      int i, rc=0;
  865.      
  866.      float ticks[N__TICKS];
  867.      char labels[N__TICKS][10];
  868.      float magnitude,y,yr,startTick,tickSize;
  869.      float x=0.0;
  870.      char string[70];
  871.      float fontSize;
  872.      float pmag;
  873.      char pstr[10];
  874.      int nLogTicks;
  875.      float logTicks[5];
  876.      float maglow, maghigh, magrng, magStep;
  877.      
  878.      axis_t *thisaxis;
  879.      char xalign='c', yalign='c';
  880.      
  881.      fontSize = MIN(disp->drawRect.size.width*0.04,14.0);
  882.  
  883.      switch (axis)
  884.      {
  885.      case XAXIS:
  886.       thisaxis = &disp->xAxis;
  887.       break;
  888.      case YAXIS:
  889.       thisaxis = &disp->yAxis;
  890.       break;
  891.      default:
  892.       return -1;
  893.      }
  894.      
  895.      /* draw ticks on y axis */
  896.      if (thisaxis->low >= thisaxis->high)
  897.      {
  898.       h_error("h_plot: axis low > high");
  899.       sprintf(string,"low = %f, high = %f\n",thisaxis->low,
  900.           thisaxis->high);
  901.       h_error(string);
  902.       return -1;
  903.      }
  904.  
  905.      if (!thisaxis->flags.log)
  906.      {
  907.       tickSize = calcTicks((thisaxis->high -
  908.                 thisaxis->low), &magnitude);
  909.       startTick = ceil( thisaxis->low / tickSize) * tickSize;
  910.  
  911.       if (fabs(magnitude) <= 3.0)
  912.            pmag = 0.0;
  913.       else 
  914.       {
  915.            if (startTick != 0.0)
  916.             pmag = floor(log10(fabs(startTick)));
  917.            else
  918.             pmag = magnitude;
  919.       }
  920.                        
  921.       sprintf(pstr,"%%1.%df",(int)MAX((pmag-magnitude),0.0));
  922.  
  923.       y = startTick;
  924.       while (y <= thisaxis->high*(1.0+NUM_FUZZ))
  925.       {
  926.            if (nticks >= N__TICKS)
  927.            {
  928.             h_error("Too many ticks along y axis.");
  929.             return -1;
  930.            }
  931.            
  932.            yr = floor(y/pow(10.0,magnitude) + 0.5);
  933.            ticks[nticks] = yr * pow(10.0,magnitude);
  934.            sprintf(labels[nticks],pstr,yr*pow(10.0,magnitude-pmag));
  935.            nticks++;
  936.            y += tickSize;
  937.       }
  938.       if (fabs(magnitude) <= 3.0) magnitude = 0.0;
  939.      }
  940.      else
  941.      {
  942.       maghigh = ceil(log10(thisaxis->high));
  943.       maglow = floor(log10(thisaxis->low));
  944.       magrng = maghigh - maglow;
  945.       
  946.       if (magrng <=3)
  947.       {
  948.            nLogTicks = 3;
  949.            logTicks[0] = 1.0;
  950.            logTicks[1] = 2.0;
  951.            logTicks[2] = 5.0;
  952.            magStep = 1.0;
  953.       }
  954.       else 
  955.       {
  956.            nLogTicks = 1;
  957.            logTicks[0] = 1.0;
  958.            if (magrng <= 7) 
  959.             magStep = 1.0;
  960.            else 
  961.             magStep = 2.0;
  962.       }
  963.            
  964.       if (nLogTicks == 3 && (fabs(maglow)>3 || fabs(maghigh)>3))
  965.            pmag = maglow;
  966.       else
  967.            pmag = 0.0;
  968.       magnitude = maglow;
  969.       i = 0;
  970.       while ((y=logTicks[i]*pow(10.0,magnitude)) <
  971.          thisaxis->high*(1.0+NUM_FUZZ))
  972.       {
  973.            if (nticks >= N__TICKS)
  974.            {
  975.             h_error("Too many ticks along y axis.");
  976.             return -1;
  977.            }
  978.            if (y >= thisaxis->low)
  979.            {
  980.             ticks[nticks] = log10(y);
  981.             /*
  982.              * be careful: there is a bug in the NeXT (s)printf 
  983.              *   routine when you do, eg. printf("%1.0g",0.01);
  984.              */
  985.             if ((magnitude-pmag) > 4 || (magnitude-pmag) < -3)
  986.              strcpy(pstr,"%1.0e");
  987.             else
  988.              sprintf(pstr,"%%1.%df",
  989.                  (int)((magnitude-pmag)>0?
  990.                        0:-(magnitude-pmag)));
  991.             sprintf(labels[nticks],pstr,y*pow(10.0,-pmag));
  992.             nticks++;
  993.            }
  994.            
  995.            i++;
  996.            if (i>=nLogTicks)
  997.            {
  998.             i = 0;
  999.             magnitude += magStep;
  1000.            }
  1001.       }
  1002.      }
  1003.      
  1004.      /*
  1005.       * Plot the tick marks
  1006.       */
  1007.      switch (axis)
  1008.      {
  1009.      case XAXIS:
  1010.       drawXTicks(ticks,nticks,thisaxis->tickLength,
  1011.              thisaxis->flags.tickLocation);
  1012.       break;
  1013.      case YAXIS:
  1014.       drawYTicks(ticks,nticks,thisaxis->tickLength,
  1015.              thisaxis->flags.tickLocation);
  1016.       break;
  1017.      default:
  1018.       return -1;
  1019.      }
  1020.      if (rc != 0) return rc;
  1021.  
  1022.      switch (axis)
  1023.      {
  1024.      case XAXIS:
  1025.       y = disp->marginRect.origin.y - YPADDING(disp);
  1026.       fontSize = MIN(fontSize, 
  1027.              (disp->marginRect.origin.y -
  1028.               disp->drawRect.origin.y)/2.0);
  1029.       xalign = 'c';
  1030.       yalign = 't';
  1031.       break;
  1032.      case YAXIS:
  1033.       x = disp->marginRect.origin.x - XPADDING(disp);
  1034.       yr = MIN(disp->drawRect.size.width*0.04,14.0);
  1035.       yr = MIN(fontSize, 
  1036.             disp->drawRect.size.height/strlen(disp->yAxis.label)*2);
  1037.       fontSize = MIN(fontSize, 
  1038.              (disp->marginRect.origin.x -
  1039.               disp->drawRect.origin.x - yr)
  1040.              /strlen(labels[nticks-1]));
  1041.       xalign = 'r';
  1042.       yalign = 'c';
  1043.       break;
  1044.      default:
  1045.       return -1;
  1046.      }
  1047.  
  1048.      for (i=0; i<nticks; i++ )
  1049.      {
  1050.       if (!thisaxis->flags.log)
  1051.       {
  1052.            switch (axis)
  1053.            {
  1054.            case XAXIS:
  1055.             x = disp->marginRect.origin.x +
  1056.              (ticks[i]-thisaxis->low) * 
  1057.                   disp->marginRect.size.width
  1058.                    /(thisaxis->high - thisaxis->low);
  1059.             break;
  1060.            case YAXIS:
  1061.             y = disp->marginRect.origin.y +
  1062.              (ticks[i]-thisaxis->low) * 
  1063.                   disp->marginRect.size.height 
  1064.                    /(thisaxis->high - thisaxis->low);
  1065.             break;
  1066.            default:
  1067.             return -1;
  1068.            }
  1069.       }
  1070.       else
  1071.       {
  1072.            switch (axis)
  1073.            {
  1074.            case XAXIS:
  1075.             x = disp->marginRect.origin.x +
  1076.              (ticks[i]-log10(thisaxis->low))
  1077.                   * disp->marginRect.size.width
  1078.                    /log10(thisaxis->high/thisaxis->low);
  1079.             break;
  1080.            case YAXIS:
  1081.             y = disp->marginRect.origin.y +
  1082.              (ticks[i]-log10(thisaxis->low))
  1083.                   * disp->marginRect.size.height 
  1084.                    /log10(thisaxis->high/thisaxis->low);
  1085.             break;
  1086.            default:
  1087.             return -1;
  1088.            }
  1089.       }
  1090.  
  1091.       if (drawText(labels[i], x, y, fontSize, 0.0,
  1092.                xalign, yalign ) != 0)
  1093.       {
  1094.            h_error("Error draw label on axis.");
  1095.            return -1;
  1096.       }
  1097.      }      
  1098.  
  1099.      if (pmag != 0.0)
  1100.      {
  1101.       switch (axis)
  1102.       {
  1103.       case XAXIS:
  1104.            x = disp->marginRect.origin.x + disp->marginRect.size.width
  1105.             + XPADDING(disp);
  1106.            y = disp->marginRect.origin.y;
  1107.            break;
  1108.       case YAXIS:
  1109.            x = disp->marginRect.origin.x;
  1110.            y = disp->marginRect.origin.y + disp->marginRect.size.height
  1111.             + YPADDING(disp);
  1112.            break;
  1113.       default:
  1114.            return -1;
  1115.       }
  1116.       drawMag(x, y, (int)pmag, fontSize*0.75 );
  1117.      }
  1118.      
  1119.      return 0;
  1120. }
  1121.  
  1122. static int drawXTicks(float *ticks, int nticks, float tickLength, char tickLoc )
  1123. {
  1124.      int rc = 0;
  1125.      
  1126.      if (tickLoc & PLOTBOTTOM)
  1127.      {
  1128.       PlotSwitch(plot_drvr,drawXTicks,(ticks,nticks,tickLength,0));
  1129.      }
  1130.      if (tickLoc & PLOTTOP)
  1131.      {
  1132.       PlotSwitch(plot_drvr,drawXTicks,(ticks,nticks,tickLength,1));
  1133.      }
  1134.      return rc;
  1135. }
  1136.  
  1137.  
  1138. static int drawYTicks(float *ticks, int nticks, float tickLength, char tickLoc )
  1139. {
  1140.      int rc = 0;
  1141.      
  1142.      if (tickLoc & PLOTLEFT)
  1143.      {
  1144.       PlotSwitch(plot_drvr,drawYTicks,(ticks,nticks,tickLength,0));
  1145.      }
  1146.      if (tickLoc & PLOTRIGHT)
  1147.      {
  1148.       PlotSwitch(plot_drvr,drawYTicks,(ticks,nticks,tickLength,1));
  1149.      }
  1150.      return rc;
  1151. }
  1152.  
  1153.  
  1154. static int drawMag( float x, float y, int mag, float fontSize )
  1155. {
  1156.      int rc;
  1157.  
  1158.      PlotSwitch(plot_drvr,drawMag,( x, y, mag, fontSize ));
  1159.      
  1160.      return rc;
  1161. }
  1162.      
  1163.  
  1164.  
  1165. static int drawLego1D(display disp, float *xylego, int *over, 
  1166.               int npts)
  1167. {
  1168.      int rc = 0;
  1169.      int end;
  1170.      int i_over;
  1171.      float ylow, yhigh;
  1172.      
  1173.      if (disp->yAxis.flags.log)
  1174.      {
  1175.       yhigh = log10(disp->yAxis.high);
  1176.       ylow = log10(disp->yAxis.low);
  1177.      }
  1178.      else
  1179.      {
  1180.       yhigh = disp->yAxis.high;
  1181.       ylow = disp->yAxis.low;
  1182.      }
  1183.      
  1184.      /*
  1185.       * first, check the first group of points, since
  1186.       * it is not associated with a bin. (last group is check
  1187.       * by last element in over)
  1188.       */
  1189.      if (xylego[1] > yhigh) xylego[1] = yhigh;
  1190.      if (xylego[1] < ylow) xylego[1] = ylow;
  1191.      if (xylego[3] > yhigh) xylego[3] = yhigh;
  1192.      if (xylego[3] < ylow) xylego[3] = ylow;
  1193.           
  1194.      i_over = 0;
  1195.      do
  1196.      {
  1197.       end = over[i_over++];
  1198.       if (xylego[4*(end+1)+1] > yhigh) xylego[4*(end+1)+1] = yhigh;
  1199.       if (xylego[4*(end+1)+1] < ylow) xylego[4*(end+1)+1] = ylow;
  1200.       if (xylego[4*(end+1)+3] > yhigh) xylego[4*(end+1)+3] = yhigh;
  1201.       if (xylego[4*(end+1)+3] < ylow) xylego[4*(end+1)+3] = ylow;
  1202.      }
  1203.      while (end < npts);
  1204.  
  1205.  
  1206.      PlotSwitch(plot_drvr,drawLine,( xylego, 2*(npts+2), SOLID ));
  1207.  
  1208.      return rc;
  1209. }
  1210.  
  1211.  
  1212. static int drawPoint1D(display disp, float *xy, int *over, int npts)
  1213. {
  1214.      int rc=0;
  1215.      int start = -1, end;
  1216.      int i_over=0;
  1217.           
  1218.      float xlow, ylow;
  1219.      float xhigh, yhigh;
  1220.      
  1221.      if (disp->xAxis.flags.log)
  1222.      {
  1223.       xhigh = log10(disp->xAxis.high);
  1224.       xlow = log10(disp->xAxis.low);
  1225.      }
  1226.      else
  1227.      {
  1228.       xhigh = disp->xAxis.high;
  1229.       xlow = disp->xAxis.low;
  1230.      }
  1231.      if (disp->yAxis.flags.log)
  1232.      {
  1233.       yhigh = log10(disp->yAxis.high);
  1234.       ylow = log10(disp->yAxis.low);
  1235.      }
  1236.      else
  1237.      {
  1238.       yhigh = disp->yAxis.high;
  1239.       ylow = disp->yAxis.low;
  1240.      }
  1241.      
  1242.      do
  1243.      {
  1244.       /*
  1245.        * Check that points marked as over are really over (it may be
  1246.        * that error bar is over).
  1247.        */
  1248.       do 
  1249.       {
  1250.            if (start == -1)
  1251.             start = 0;
  1252.            else
  1253.             start = over[i_over++]+1;
  1254.            end = over[i_over];
  1255.  
  1256.            while (end<npts)
  1257.            {
  1258.             /* include some fuzz in comparisons */
  1259.             if (xy[2*end+1]<(ylow*(1.0-NUM_FUZZ)) || 
  1260.             xy[2*end+1]>(yhigh*(1.0+NUM_FUZZ))) break;
  1261.             if (xy[2*end]<(xlow*(1.0-NUM_FUZZ)) || 
  1262.             xy[2*end]>(xhigh*(1.0+NUM_FUZZ))) break;
  1263.             end = over[++i_over];
  1264.            }
  1265.       }
  1266.       while ( (end-start) <= 0 && start<npts );
  1267.             
  1268.       if (start >= npts) continue;
  1269.       
  1270.       PlotSwitch(plot_drvr,drawPoints,(&(xy[2*start]),end-start,
  1271.                        disp->plotSymbol,disp->symbolSize));
  1272.      }
  1273.      while (rc == 0 && end < npts && start < npts);
  1274.      
  1275.      return rc;
  1276. }
  1277.  
  1278.  
  1279. static int drawLine1D(display disp, float *xy, int *over, 
  1280.               int npts, int xmin_pt, linestyle_t ls)
  1281. {
  1282.      int rc = 0;
  1283.      int start, end;
  1284.      int s_saved, e_saved;
  1285.      int i_over = 0;
  1286.      float *xy_alt;
  1287.      int *over_alt;
  1288.      float xlow, ylow;
  1289.      float xhigh, yhigh;
  1290.      
  1291.      float start_save[2],end_save[2];
  1292.  
  1293.      if (npts < 2) return 0;
  1294.      
  1295.      if (xmin_pt > 0)
  1296.      {
  1297.       if ( (xy_alt = (float *)malloc( 2*npts*sizeof(float) ))
  1298.           == NULL)
  1299.       {
  1300.            h_error("drawLine1D - ran out of memory");
  1301.            return -1;
  1302.       }
  1303.       if ( (over_alt = (int *)malloc( (npts+1)*sizeof(int) ))
  1304.           == NULL)
  1305.       {
  1306.            h_error("drawLine1D - ran out of memory");
  1307.            return -1;
  1308.       }
  1309.       memcpy((char *)xy_alt, (char *)&xy[2*xmin_pt],
  1310.          2*(npts-xmin_pt)*sizeof(float) );
  1311.       memcpy((char *)&xy_alt[2*(npts-xmin_pt)], (char *)xy,
  1312.          2*xmin_pt*sizeof(float) );
  1313.       for (i_over=0; over[i_over]<npts; i_over++ )
  1314.       {
  1315.            if (over[i_over] >= xmin_pt)
  1316.             over_alt[i_over] = over[i_over] - xmin_pt;
  1317.            else
  1318.             over_alt[i_over] = over[i_over] + npts - xmin_pt;
  1319.       }
  1320.       over_alt[i_over++] = npts;
  1321.       if (i_over > 1)
  1322.            qsort( over_alt, (size_t)i_over, sizeof(int), intcmp );
  1323.      }
  1324.      else
  1325.      {
  1326.       over_alt = over;
  1327.       xy_alt = xy;
  1328.      }
  1329.           
  1330.      if (disp->xAxis.flags.log)
  1331.      {
  1332.       xhigh = log10(disp->xAxis.high);
  1333.       xlow = log10(disp->xAxis.low);
  1334.      }
  1335.      else
  1336.      {
  1337.       xhigh = disp->xAxis.high;
  1338.       xlow = disp->xAxis.low;
  1339.      }
  1340.      if (disp->yAxis.flags.log)
  1341.      {
  1342.       yhigh = log10(disp->yAxis.high);
  1343.       ylow = log10(disp->yAxis.low);
  1344.      }
  1345.      else
  1346.      {
  1347.       yhigh = disp->yAxis.high;
  1348.       ylow = disp->yAxis.low;
  1349.      }
  1350.      
  1351.      i_over = 0;
  1352.      start = -1;
  1353.      do
  1354.      {
  1355.       /*
  1356.        * Check that points marked as over are really over (it may be
  1357.        * that error bar is over). Interpolate to find point where
  1358.        * line crosses axes.
  1359.        */
  1360.       s_saved = 0;
  1361.       e_saved = 0;
  1362.  
  1363.       do 
  1364.       {
  1365.            if (start == -1)
  1366.             start = 0;
  1367.            else
  1368.             start = over_alt[i_over++]+1;
  1369.            end = over_alt[i_over];
  1370.  
  1371.            while (end<npts)
  1372.            {
  1373.             /* include some fuzz in comparisons */
  1374.             if (xy_alt[2*end+1]<(ylow*(1.0-NUM_FUZZ)) || 
  1375.             xy_alt[2*end+1]>(yhigh*(1.0+NUM_FUZZ))) break;
  1376.             if (xy_alt[2*end]<(xlow*(1.0-NUM_FUZZ)) ||
  1377.             xy_alt[2*end]>(xhigh*(1.0+NUM_FUZZ))) break;
  1378.             end = over_alt[++i_over];
  1379.            }
  1380.       }
  1381.       while ( (end-start) <= 0 && start<npts );
  1382.             
  1383.       if (start >= npts) continue;
  1384.  
  1385.       if (start > 0)
  1386.       {
  1387.            memcpy(start_save,&(xy_alt[2*(start-1)]),2*sizeof(float) );
  1388.            s_saved = 1;
  1389.            box_cross( &xy_alt[2*(start-1)], &xy_alt[2*start], disp);
  1390.            start--;
  1391.       }
  1392.       if (end < npts)
  1393.       {
  1394.            memcpy(end_save,&(xy_alt[2*end]),2*sizeof(float) );
  1395.            e_saved = 1;
  1396.            box_cross( &xy_alt[2*end], &xy_alt[2*(end-1)], disp);
  1397.            end++;
  1398.       }
  1399.  
  1400.       PlotSwitch(plot_drvr,drawLine,
  1401.              (&(xy_alt[2*start]),end-start,ls));
  1402.  
  1403.       if (s_saved)
  1404.            memcpy(&(xy_alt[2*start]),start_save,2*sizeof(float) );
  1405.       if (e_saved)
  1406.            memcpy(&(xy_alt[2*(end-1)]),end_save,2*sizeof(float) );
  1407.  
  1408.  
  1409.      }
  1410.      while (rc == 0 && end < npts && start < npts);
  1411.  
  1412.      if (xmin_pt != 0)
  1413.      {
  1414.       free(xy_alt);
  1415.       free(over_alt);
  1416.      }
  1417.      
  1418.      return rc;
  1419. }
  1420.  
  1421. /*
  1422.  * integer comparison for qsort.
  1423.  */
  1424. int intcmp(const void *i1, const void *i2)
  1425. {
  1426.      if ((int *)i1 > (int *)i2)
  1427.       return 1;
  1428.      else if ((int *)i1 == (int *)i2)
  1429.       return 0;
  1430.      else 
  1431.       return -1;
  1432. }
  1433.  
  1434. /*
  1435.  * find point where interpolation between p1 and p2 cross the 
  1436.  *  axes box.
  1437.  */
  1438. static void box_cross( float *p1, float *p2, display disp )
  1439. {
  1440.      float xlow,ylow, xhigh, yhigh;
  1441.      float x_int1,x_int2,y_int1,y_int2;
  1442.      
  1443.      if (disp->xAxis.flags.log)
  1444.      {
  1445.       xlow = log10(disp->xAxis.low);
  1446.       xhigh = log10(disp->xAxis.high);
  1447.      }
  1448.      else
  1449.      {
  1450.       xlow = disp->xAxis.low;
  1451.       xhigh = disp->xAxis.high;
  1452.      }
  1453.      if (disp->yAxis.flags.log)
  1454.      {
  1455.       ylow = log10(disp->yAxis.low);
  1456.       yhigh = log10(disp->yAxis.high);
  1457.      }
  1458.      else
  1459.      {
  1460.       ylow = disp->yAxis.low;
  1461.       yhigh = disp->yAxis.high;
  1462.      }
  1463.      
  1464.      x_int1 = p1[0] + (ylow-p1[1])*(p2[0]-p1[0])/(p2[1]-p1[1]);
  1465.      if (x_int1 >= xlow && x_int1 <= xhigh &&
  1466.      ylow >= p1[1] && ylow <= p2[1])
  1467.      {
  1468.       p1[0] = x_int1;
  1469.       p1[1] = ylow;
  1470.       return;
  1471.      }
  1472.  
  1473.      x_int2 = p1[0] + (yhigh-p1[1])*(p2[0]-p1[0])/(p2[1]-p1[1]);
  1474.      if (x_int2 >= xlow && x_int2 <= xhigh &&
  1475.      yhigh <= p1[1] && yhigh >= p2[1])
  1476.      {
  1477.       p1[0] = x_int2;
  1478.       p1[1] = yhigh;
  1479.       return;
  1480.      }
  1481.  
  1482.      y_int1 = p1[1] + (xlow-p1[0])*(p2[1]-p1[1])/(p2[0]-p1[0]);
  1483.      if (y_int1 >= ylow && y_int1 <= yhigh &&
  1484.      xlow >= p1[0] && xlow <= p2[0])
  1485.      {
  1486.       p1[0] = xlow;
  1487.       p1[1] = y_int1;
  1488.       return;
  1489.      }
  1490.  
  1491.      y_int2 = p1[1] + (xhigh-p1[0])*(p2[1]-p1[1])/(p2[0]-p1[0]);
  1492.      if (y_int2 >= ylow && y_int1 <= yhigh &&
  1493.      xhigh <= p1[0] && xhigh >= p2[0])
  1494.      {
  1495.       p1[0] = xhigh;
  1496.       p1[1] = y_int2;
  1497.       return;
  1498.      }
  1499. }
  1500.  
  1501.  
  1502. /*
  1503.  * draw errors - err_type specifies whether they are x or y errors.
  1504.  */
  1505. static int drawError1D(display disp, float *xy, float *err, 
  1506.                int *over, int npts, int err_type)
  1507. {
  1508.      int rc = 0;
  1509.      int start = -1;
  1510.      int end;
  1511.      int i_over = 0;
  1512.      
  1513.      float xlow, ylow;
  1514.      float xhigh, yhigh;
  1515.      
  1516.      if (disp->xAxis.flags.log)
  1517.      {
  1518.       xhigh = log10(disp->xAxis.high);
  1519.       xlow = log10(disp->xAxis.low);
  1520.      }
  1521.      else
  1522.      {
  1523.       xhigh = disp->xAxis.high;
  1524.       xlow = disp->xAxis.low;
  1525.      }
  1526.      if (disp->yAxis.flags.log)
  1527.      {
  1528.       yhigh = log10(disp->yAxis.high);
  1529.       ylow = log10(disp->yAxis.low);
  1530.      }
  1531.      else
  1532.      {
  1533.       yhigh = disp->yAxis.high;
  1534.       ylow = disp->yAxis.low;
  1535.      }
  1536.      
  1537.      do
  1538.      {
  1539.       /*
  1540.        * Plot any x/y error bar where the y/x coordinate is in range
  1541.        *  because of bars which reach into plot, even when point
  1542.        *  is outside.
  1543.        */
  1544.       do 
  1545.       {
  1546.            if (start == -1)
  1547.             start = 0;
  1548.            else
  1549.             start = over[i_over++]+1;
  1550.            end = over[i_over];
  1551.  
  1552.            while (end<npts)
  1553.            {
  1554.             if (err_type==XERROR)
  1555.             {
  1556.              /* include some fuzz in comparisons */
  1557.              if ( xy[2*end+1]<(ylow*(1.0-NUM_FUZZ)) || 
  1558.                  xy[2*end+1]>(yhigh*(1.0+NUM_FUZZ))) break;
  1559.              if ( err[2*end]<(xlow*(1.0-NUM_FUZZ)) || 
  1560.                  err[2*end+1]>(xhigh*(1.0+NUM_FUZZ))) break;
  1561.              if ( err[2*end] > xhigh)
  1562.                   err[2*end] = xhigh;
  1563.              if ( err[2*end+1] < xlow)
  1564.                   err[2*end+1] = xlow;
  1565.             }
  1566.             else
  1567.             {
  1568.              if ( xy[2*end]<(xlow*(1.0-NUM_FUZZ)) || 
  1569.                  xy[2*end]>(xhigh*(1.0+NUM_FUZZ))) break;
  1570.              if ( err[2*end]<(ylow*(1.0-NUM_FUZZ)) ||
  1571.                  err[2*end+1]>(yhigh*(1.0+NUM_FUZZ))) break;
  1572.              if ( err[2*end] > yhigh)
  1573.                   err[2*end] = yhigh;
  1574.              if ( err[2*end+1] < ylow)
  1575.                   err[2*end+1] = ylow;
  1576.             }
  1577.             end = over[++i_over];
  1578.            }
  1579.       }
  1580.       while ( (end-start) <= 0 && start<npts );
  1581.             
  1582.       if (start >= npts) continue;
  1583.       
  1584.       if (err_type == YERROR)
  1585.       {
  1586.            PlotSwitch(plot_drvr,drawYError,(&(xy[2*start]),
  1587.                         &(err[2*start]),
  1588.                         end-start));
  1589.       }
  1590.       else
  1591.       {
  1592.            PlotSwitch(plot_drvr,drawXError,(&(xy[2*start]),
  1593.                         &(err[2*start]),
  1594.                         end-start));
  1595.       }
  1596.      }
  1597.      while (rc == 0 && end < npts && start < npts);
  1598.      
  1599.      return rc;
  1600. }
  1601.  
  1602.  
  1603. static int drawColor2D(display disp)
  1604. {
  1605.      int rc=0;
  1606.  
  1607.      if (!init_done)
  1608.      {
  1609.       initPlot( disp );
  1610.       init_done = 1;
  1611.      }
  1612.           
  1613.      rc = disp->drawtype & COLOR;      /* use temporarily */
  1614.      switch (plot_drvr)
  1615.      {
  1616.       NeXTFunc(drawColor2D,(disp->bins.xAxis.nBins,
  1617.                 disp->bins.yAxis.nBins,
  1618.                 &disp->marginRect,
  1619.                 disp->bins.data,
  1620.                 disp->bins.binMin,
  1621.                 disp->bins.binMax,
  1622.                 rc));
  1623.  
  1624.       PSFunc(drawColor2D,(disp->bins.xAxis.nBins,
  1625.                   disp->bins.yAxis.nBins,
  1626.                   &disp->marginRect,
  1627.                   disp->bins.data,
  1628.                   disp->bins.binMin,
  1629.                   disp->bins.binMax,
  1630.                   rc));
  1631.  
  1632.       UPFunc(drawColor2D,(disp));
  1633.       
  1634.       XIVFunc(drawColor2D,(disp->bins.xAxis.nBins,
  1635.                    disp->bins.yAxis.nBins,
  1636.                    disp->bins.binMin,
  1637.                    disp->bins.binMax,
  1638.                    disp->bins.data));
  1639.  
  1640.       X11Func(drawColor2D,(disp->bins.xAxis.nBins,
  1641.                    disp->bins.yAxis.nBins,
  1642.                    disp->bins.binMin,
  1643.                    disp->bins.binMax,
  1644.                    disp->bins.data,
  1645.                    rc));
  1646.  
  1647.      default:
  1648.       rc = -1;
  1649.      }
  1650.      
  1651.      return rc;
  1652. }
  1653.  
  1654.  
  1655. static int drawScatter2D(display disp)
  1656. {
  1657.      int rc;
  1658.      int npts;
  1659.      float *xylist;
  1660.      struct maxs max;
  1661.      
  1662.      npts = getScatData( disp, &xylist, &max );
  1663.      if (npts == 0) return 0;
  1664.      if (npts < 0) 
  1665.      {
  1666.       h_error("h_plot - error getting scatter plot data");
  1667.       return -1;
  1668.      }
  1669.  
  1670.      if (disp->xAxis.flags.autoScale)
  1671.      {
  1672.       disp->xAxis.low = max.xl;
  1673.       disp->xAxis.high = max.xh;
  1674.       h_adjustAxis(&(disp->xAxis.low),&(disp->xAxis.high),
  1675.                0, disp->xAxis.flags.log);
  1676.      }
  1677.      if (disp->yAxis.flags.autoScale)
  1678.      {
  1679.       disp->yAxis.low = max.yl;
  1680.       disp->yAxis.high = max.yh;
  1681.       h_adjustAxis(&(disp->yAxis.low),&(disp->yAxis.high),
  1682.                0, disp->yAxis.flags.log);
  1683.      }
  1684.      
  1685.      if (!init_done)
  1686.      {
  1687.       initPlot( disp );
  1688.       init_done = 1;
  1689.      }
  1690.           
  1691.      PlotSwitch(plot_drvr,drawPoints,(xylist,npts,disp->plotSymbol,
  1692.                       disp->symbolSize));
  1693.      
  1694.      free( xylist );
  1695.      return rc;
  1696. }
  1697.  
  1698.  
  1699. static int getXYData( display disp, float **xy, float **xerr, 
  1700.              float **yerr, int **over, struct maxs *m, int *xmin_pt)
  1701. {
  1702.      int i;
  1703.      float *xp,*yp,*nt;
  1704.      float *last;
  1705.      float xlow, xhigh;
  1706.      float ylow, yhigh;
  1707.      float *xerrorp = NULL, *yerrorp = NULL;
  1708.      float *xerrp, *yerrp;
  1709.      float *xyp;
  1710.      int *overp, isover;
  1711.      int nt_dim = disp->ntuple->ndim;
  1712.      int ylog, xlog;
  1713.      
  1714.  
  1715.      if (disp->binding.x < 0 || disp->binding.y < 0 ||
  1716.      disp->binding.x >= disp->ntuple->ndim ||
  1717.      disp->binding.y >= disp->ntuple->ndim )
  1718.       return -1;
  1719.      
  1720.      m->xl = FLT_MAX;
  1721.      m->xh = -FLT_MAX;
  1722.      m->yl = FLT_MAX;
  1723.      m->yh = -FLT_MAX;
  1724.      
  1725.      ylog = disp->yAxis.flags.log;
  1726.      xlog = disp->xAxis.flags.log;
  1727.      
  1728.      if (disp->xAxis.flags.autoScale)
  1729.      {
  1730.       xlow = -FLT_MAX;
  1731.       xhigh = FLT_MAX;
  1732.      }
  1733.      else
  1734.      {
  1735.       xlow = disp->xAxis.low;
  1736.       xhigh = disp->xAxis.high;
  1737.      }
  1738.      if (disp->yAxis.flags.autoScale)
  1739.      {
  1740.       ylow = -FLT_MAX;
  1741.       yhigh = FLT_MAX;
  1742.      }
  1743.      else
  1744.      {
  1745.       ylow = disp->yAxis.low;
  1746.       yhigh = disp->yAxis.high;
  1747.      }
  1748.      
  1749.      
  1750.      memset((char *)disp->bins.totals, (char)0, sizeof(disp->bins.totals) );
  1751.  
  1752.      if ( (*xy = (float *)malloc( 2*disp->ntuple->ndata*sizeof(float)))
  1753.      == NULL )
  1754.      {
  1755.       h_error("h_plot has run out of memory.");
  1756.       return -1;
  1757.      }
  1758.      if ( disp->drawtype&ERRBAR && disp->binding.xerror >=0 )
  1759.      {
  1760.       if ((*xerr = (float *)malloc( 2*disp->ntuple->ndata*sizeof(float)))
  1761.           == NULL )
  1762.       {
  1763.            free(*xy);
  1764.            h_error("h_plot has run out of memory.");
  1765.            return -1;
  1766.       }
  1767.      }
  1768.      else *xerr = NULL;
  1769.      
  1770.      if ( disp->drawtype&ERRBAR && disp->binding.yerror >= 0 )
  1771.      {
  1772.       if ( (*yerr = (float *)malloc( 2*disp->ntuple->ndata*sizeof(float)))
  1773.           == NULL )
  1774.       {
  1775.            free(xy);
  1776.            free(xerr);
  1777.            h_error("h_plot has run out of memory.");
  1778.            return -1;
  1779.       }
  1780.      }
  1781.      else *yerr = NULL;
  1782.      
  1783.      if ( (*over = (int *)malloc( (disp->ntuple->ndata+1)*sizeof(int)))
  1784.      == NULL )
  1785.      {
  1786.       free(xy);
  1787.       free(xerr);
  1788.       free(yerr);
  1789.       h_error("h_plot has run out of memory.");
  1790.       return -1;
  1791.      }
  1792.      
  1793.      last = &( disp->ntuple->data[INDEX(disp->ntuple->ndata,
  1794.                     0,
  1795.                     disp->ntuple->ndim)] );
  1796.      xp = &(disp->ntuple->data[INDEX(0,disp->binding.x,disp->ntuple->ndim)]);
  1797.      yp = &(disp->ntuple->data[INDEX(0,disp->binding.y,disp->ntuple->ndim)]);
  1798.      if (disp->binding.xerror >=0 )
  1799.       xerrorp = &(disp->ntuple->data[INDEX(0,disp->binding.xerror,
  1800.                            disp->ntuple->ndim)]);
  1801.      if (disp->binding.yerror >=0 )
  1802.       yerrorp = &(disp->ntuple->data[INDEX(0,disp->binding.yerror,
  1803.                            disp->ntuple->ndim)]);
  1804.  
  1805.      nt = &(disp->ntuple->data[INDEX(0,0,disp->ntuple->ndim)]);
  1806.      
  1807.      i = 0;
  1808.      xerrp = *xerr;
  1809.      yerrp = *yerr;
  1810.      xyp = *xy;
  1811.      overp = *over;
  1812.      
  1813.      for ( ; xp<last; 
  1814.       xp += nt_dim, yp += nt_dim, nt += nt_dim,
  1815.       xerrorp += nt_dim, yerrorp += nt_dim )
  1816.      {
  1817.       if (disp->nt_cut!=NULL && doCuts(nt, disp->nt_cut)) continue;
  1818.       
  1819.       isover = 0;
  1820.       
  1821.       if (*xp > xhigh)
  1822.       {
  1823.            if (*yp > yhigh)
  1824.             disp->bins.totals[2][2] += 1.0;
  1825.            else if (*yp < ylow)
  1826.             disp->bins.totals[2][0] += 1.0;
  1827.            else
  1828.             disp->bins.totals[2][1] += 1.0;
  1829.       }
  1830.       else if (*xp < xlow)
  1831.       {
  1832.            if (*yp > yhigh)
  1833.             disp->bins.totals[0][2] += 1.0;
  1834.            else if (*yp < ylow)
  1835.             disp->bins.totals[0][0] += 1.0;
  1836.            else
  1837.             disp->bins.totals[0][1] += 1.0;
  1838.       }
  1839.       else
  1840.       {
  1841.            if (*yp > yhigh)
  1842.             disp->bins.totals[1][2] += 1.0;
  1843.            else if (*yp < ylow)
  1844.             disp->bins.totals[1][0] += 1.0;
  1845.            else
  1846.             disp->bins.totals[1][1] += 1.0;
  1847.       }
  1848.       
  1849.       if (!xlog)
  1850.            *xyp++ = *xp;
  1851.       else
  1852.       {
  1853.            if (*xp > 0.0)
  1854.             *xyp++ = log10(*xp);
  1855.            else
  1856.            {
  1857.             isover = 1;
  1858.             *xyp++ = -FLT_MAX;
  1859.            }
  1860.       }
  1861.       if (*xp > xhigh || *xp < xlow) isover = 1;
  1862.       
  1863.       if (!ylog)
  1864.            *xyp++ = *yp;
  1865.       else
  1866.       {
  1867.            if (*yp > 0.0)
  1868.             *xyp++ = log10(*yp);
  1869.            else
  1870.            {
  1871.             isover = 1;
  1872.             *xyp++ = -FLT_MAX;
  1873.            }
  1874.       }
  1875.       if (*yp > yhigh || *yp < ylow ) isover = 1;
  1876.  
  1877.       if (! isover)
  1878.       {
  1879.            if (*xp > m->xh) m->xh = *xp;
  1880.            if (*xp < m->xl) 
  1881.            {
  1882.             m->xl = *xp;
  1883.             *xmin_pt = i;
  1884.            }
  1885.            
  1886.            if (*yp > m->yh) m->yh = *yp;
  1887.            if (*yp < m->yl) m->yl = *yp;
  1888.       }
  1889.       
  1890.       
  1891.       if (xerrp != NULL) 
  1892.       {
  1893.            *xerrp = *xp + *xerrorp;
  1894.            if (!isover && *xerrp > m->xh) m->xh = *xerrp;
  1895.            if (xlog)
  1896.            {
  1897.             if (*xerrp > 0.0)
  1898.              *xerrp = log10(*xerrp);
  1899.             else 
  1900.             {
  1901.              *xerrp = -FLT_MAX;
  1902.              isover = 1;
  1903.             }
  1904.            }
  1905.            if (*xerrp > xhigh || *xerrp < xlow) isover = 1;
  1906.            xerrp++;
  1907.            
  1908.            *xerrp = *xp - *xerrorp;
  1909.            if (!isover && *xerrp < m->xl &&
  1910.            (!xlog || (xlog && *xerrp>0.0)) ) m->xl = *xerrp;
  1911.            if (xlog) 
  1912.            {
  1913.             if (*xerrp > 0.0)
  1914.              *xerrp = log10(*xerrp);
  1915.             else
  1916.             {
  1917.              *xerrp = -FLT_MAX;
  1918.              isover = 1;
  1919.             }
  1920.            }
  1921.            if (*xerrp > xhigh || *xerrp < xlow) isover = 1;
  1922.            xerrp++;
  1923.       }
  1924.       if (yerrp != NULL) 
  1925.       {
  1926.            *yerrp = *yp + *yerrorp;
  1927.            if (!isover && *yerrp > m->yh) m->yh = *yerrp;
  1928.            if (ylog)
  1929.            {
  1930.             if (*yerrp > 0.0)
  1931.              *yerrp = log10(*yerrp);
  1932.             else
  1933.             {
  1934.              *yerrp = -FLT_MAX;
  1935.              isover = 1;
  1936.             }
  1937.            }
  1938.            if (*yerrp > yhigh || *yerrp < ylow ) isover = 1;
  1939.            yerrp++;
  1940.            
  1941.            *yerrp = *yp - *yerrorp;
  1942.            if (!isover && *yerrp < m->yl &&
  1943.            (!ylog || (ylog && *yerrp > 0.0)) )
  1944.             m->yl = *yerrp;
  1945.            if (ylog)
  1946.            {
  1947.             if (*yerrp > 0.0)
  1948.              *yerrp = log10(*yerrp);
  1949.             else
  1950.             {
  1951.              *yerrp = -FLT_MAX;
  1952.              isover = 1;
  1953.             }
  1954.            }
  1955.            if (*yerrp > yhigh || *yerrp < ylow ) isover = 1;
  1956.            yerrp++;
  1957.       }
  1958.       
  1959.       if (isover)
  1960.       {
  1961.            *overp++ = i;
  1962.            isover = 1;
  1963.       }
  1964.  
  1965.       i++;
  1966.      }
  1967.      
  1968.      /* indicate last item in over */
  1969.      *overp = i;
  1970.      
  1971.      return i;
  1972. }
  1973.  
  1974. static int getHistoData( display disp, float **xy, float **xylego,
  1975.             float **xerr, float **yerr, int **over,
  1976.             struct maxs *m)
  1977. {
  1978.      int i;
  1979.      float binwidth, binside;
  1980.      float *yp;
  1981.      float *last;
  1982.      float *xyp;
  1983.      int *overp;
  1984.      float xlow,xhigh;
  1985.      float ylow, yhigh;
  1986.      float *yerrorp, *xerrp, *yerrp;
  1987.      float *xylegop;
  1988.      float def_bin;
  1989.      int ylog, isover;
  1990.      
  1991.      if (disp->bins.data == NULL) return -1;
  1992.  
  1993.      ylog = disp->yAxis.flags.log;
  1994.      
  1995.      if (disp->xAxis.flags.autoScale)
  1996.      {
  1997.       xlow = -FLT_MAX;
  1998.       xhigh = FLT_MAX;
  1999.      }
  2000.      else
  2001.      {
  2002.       xlow = disp->xAxis.low;
  2003.       xhigh = disp->xAxis.high;
  2004.      }
  2005.      if (disp->yAxis.flags.autoScale)
  2006.      {
  2007.       ylow = -FLT_MAX;
  2008.       yhigh = FLT_MAX;
  2009.      }
  2010.      else
  2011.      {
  2012.       ylow = disp->yAxis.low;
  2013.       yhigh = disp->yAxis.high;
  2014.      }
  2015.      
  2016.      m->xl = FLT_MAX;
  2017.      m->xh = -FLT_MAX;
  2018.      m->yl = FLT_MAX;
  2019.      m->yh = -FLT_MAX;
  2020.      
  2021.      if ( (*xy = (float *)malloc( 2*disp->bins.xAxis.nBins*sizeof(float)))
  2022.      == NULL )
  2023.      {
  2024.       h_error("h_plot has run out of memory.");
  2025.       return -1;
  2026.      }
  2027.      if ( (*xylego = (float *)malloc( 4*(disp->bins.xAxis.nBins+2)*sizeof(float)))
  2028.      == NULL )
  2029.      {
  2030.       free(*xy);
  2031.       h_error("h_plot has run out of memory.");
  2032.       return -1;
  2033.      }
  2034.      if ( disp->drawtype&ERRBAR )
  2035.      {
  2036.       if ( (*xerr = (float *)malloc( 2*disp->bins.xAxis.nBins*sizeof(float)))
  2037.           == NULL )
  2038.       {
  2039.            free(*xy);
  2040.            free(*xylego);
  2041.            h_error("h_plot has run out of memory.");
  2042.            return -1;
  2043.       }
  2044.      }
  2045.      else *xerr = NULL;
  2046.      
  2047.      if ( disp->drawtype&ERRBAR )
  2048.      {
  2049.       if ( (*yerr = (float *)malloc( 2*disp->bins.xAxis.nBins*sizeof(float)))
  2050.           == NULL )
  2051.       {
  2052.            free(*xy);
  2053.            free(*xylego);
  2054.            free(*xerr);
  2055.            h_error("h_plot has run out of memory.");
  2056.            return -1;
  2057.       }
  2058.      }
  2059.      else *yerr = NULL;
  2060.      
  2061.      if ( (*over = (int *)malloc( (disp->bins.xAxis.nBins+1)*sizeof(int)))
  2062.      == NULL )
  2063.      {
  2064.       free(*xy);
  2065.       free(*xylego);
  2066.       free(*xerr);
  2067.       free(*yerr);
  2068.       h_error("h_plot has run out of memory.");
  2069.       return -1;
  2070.      }
  2071.      
  2072.      last = &( disp->bins.data[disp->bins.xAxis.nBins] );
  2073.      yp = disp->bins.data;
  2074.      yerrorp = disp->bins.variance;
  2075.      
  2076.      i = 0;
  2077.      xerrp = *xerr;
  2078.      yerrp = *yerr;
  2079.      xyp = *xy;
  2080.      xylegop = *xylego;
  2081.      overp = *over;
  2082.      binside = disp->xAxis.low;
  2083.      binwidth = (disp->xAxis.high - disp->xAxis.low)/disp->bins.xAxis.nBins;
  2084.  
  2085.      /*
  2086.       * connect the bins to the walls of the plot.
  2087.       * this is only useful if we allow plot axis to be bigger 
  2088.       * than the bin limits.
  2089.       */
  2090.      if (disp->yAxis.low > 0.0) def_bin = disp->yAxis.low;
  2091.      else def_bin = 0.0;
  2092.      if (ylog)
  2093.       if (def_bin > 0.0)
  2094.            def_bin = log10(def_bin);
  2095.       else
  2096.            def_bin = -FLT_MAX;
  2097.      
  2098.      *xylegop++ = disp->xAxis.low;
  2099.      *xylegop++ = def_bin;
  2100.      *xylegop++ = binside;
  2101.      *xylegop++ = def_bin;
  2102.      
  2103.      for ( ; yp<last; yp++, yerrorp++ )
  2104.      {
  2105.       isover = 0;
  2106.       
  2107.       *xyp++ = binside + binwidth/2.0;
  2108.       
  2109.       if (!ylog)
  2110.            *xyp++ = *yp;
  2111.       else
  2112.       {
  2113.            if (*yp > 0.0)
  2114.             *xyp++ = log10(*yp);
  2115.            else
  2116.            {
  2117.             isover = 1;
  2118.             *xyp++ = -FLT_MAX;
  2119.            }
  2120.       }
  2121.       if (*yp > yhigh || *yp < ylow) isover = 1;
  2122.       
  2123.       if (!isover)
  2124.       {
  2125.            if (*yp > m->yh) m->yh = *yp;
  2126.            if (*yp < m->yl && (!ylog || (ylog && *yp>0.0)) )
  2127.             m->yl = *yp;
  2128.       }
  2129.       
  2130.       
  2131.       *xylegop++ = binside;
  2132.       if (xerrp != NULL) *xerrp++ = binside;
  2133.       *xylegop = *yp;
  2134.       if (ylog)
  2135.       {
  2136.            if (*xylegop > 0.0)
  2137.             *xylegop = log10(*xylegop);
  2138.            else
  2139.             *xylegop = -FLT_MAX;
  2140.       }
  2141.       xylegop++;
  2142.       
  2143.       binside += binwidth;
  2144.       *xylegop++ = binside;
  2145.       *xylegop = *(xylegop-2);
  2146.       xylegop++;
  2147.  
  2148.       if (xerrp != NULL) *xerrp++ = binside;
  2149.       
  2150.       if (yerrp != NULL) 
  2151.       {
  2152.            *yerrp = *yp + sqrt(*yerrorp);
  2153.            if (!isover && *yerrp > m->yh) m->yh = *yerrp;
  2154.            if (ylog) 
  2155.            {
  2156.             if (*yerrp > 0.0)
  2157.              *yerrp = log10(*yerrp);
  2158.             else
  2159.             {
  2160.              *yerrp = -FLT_MAX;
  2161.              isover = 1;
  2162.             }
  2163.            }
  2164.            if (*yerrp > yhigh || *yerrp < ylow) isover = 1;
  2165.            yerrp++;
  2166.  
  2167.            *yerrp = *yp - sqrt(*yerrorp);
  2168.            if (!isover && *yerrp < m->yl && 
  2169.            !( ylog && *yerrp<=0.0))
  2170.             m->yl = *yerrp;
  2171.            if (ylog)
  2172.            {
  2173.             if (*yerrp > 0.0)
  2174.              *yerrp = log10(*yerrp);
  2175.             else
  2176.             {
  2177.              *yerrp = -FLT_MAX;
  2178.              isover = 1;
  2179.             }
  2180.            }
  2181.            if (*yerrp > yhigh || *yerrp < ylow) isover = 1;
  2182.            yerrp++;
  2183.       }
  2184.  
  2185.       if (isover)
  2186.       {
  2187.            isover = 1;
  2188.            *overp++ = i;
  2189.       }
  2190.       
  2191.       i++;
  2192.      }
  2193.  
  2194.      /* connect to the walls */
  2195.      *xylegop++ = binside;
  2196.      *xylegop++ = def_bin;
  2197.      *xylegop++ = disp->xAxis.high;
  2198.      *xylegop++ = def_bin;
  2199.      
  2200.      /* indicate last item in over. */
  2201.      *overp = i;
  2202.      
  2203.      return i;
  2204. }
  2205.  
  2206.  
  2207. static int getScatData( display disp, float **xy, struct maxs *m )
  2208. {
  2209.      int i = 0;
  2210.      float *xp,*yp,*nt;
  2211.      float *last;
  2212.      float *xyp;
  2213.      int nt_dim = disp->ntuple->ndim;
  2214.      float xlow, xhigh;
  2215.      float ylow, yhigh;
  2216.      int xlog, ylog;
  2217.      
  2218.  
  2219.      if (disp->binding.x < 0 || disp->binding.y < 0 ||
  2220.      disp->binding.x >= disp->ntuple->ndim ||
  2221.      disp->binding.y >= disp->ntuple->ndim )
  2222.       return -1;
  2223.  
  2224.      xlog = disp->xAxis.flags.log;
  2225.      ylog = disp->yAxis.flags.log;
  2226.      
  2227.      m->xl = FLT_MAX;
  2228.      m->xh = -FLT_MAX;
  2229.      m->yl = FLT_MAX;
  2230.      m->yh = -FLT_MAX;
  2231.      
  2232.      if (disp->xAxis.flags.autoScale)
  2233.      {
  2234.       xlow = -FLT_MAX;
  2235.       xhigh = FLT_MAX;
  2236.      }
  2237.      else
  2238.      {
  2239.       xlow = disp->xAxis.low;
  2240.       xhigh = disp->xAxis.high;
  2241.      }
  2242.      if (disp->yAxis.flags.autoScale)
  2243.      {
  2244.       ylow = -FLT_MAX;
  2245.       yhigh = FLT_MAX;
  2246.      }
  2247.      else
  2248.      {
  2249.       ylow = disp->yAxis.low;
  2250.       yhigh = disp->yAxis.high;
  2251.      }
  2252.      
  2253.      if (xlog && xlow < FLT_MIN) xlow = FLT_MIN;
  2254.      if (ylog && ylow < FLT_MIN) ylow = FLT_MIN;
  2255.      if (xhigh <= xlow) xhigh = xlow + 1.0;
  2256.      if (yhigh <= ylow) yhigh = ylow + 1.0;
  2257.      
  2258.      if ( (*xy = (float *)malloc( 2*disp->ntuple->ndata*sizeof(float)))
  2259.      == NULL )
  2260.      {
  2261.       h_error("h_plot has run out of memory.");
  2262.       return -1;
  2263.      }
  2264.      
  2265.      memset((char *)disp->bins.totals, (char)0, sizeof(disp->bins.totals) );
  2266.  
  2267.      last = &( disp->ntuple->data[INDEX(disp->ntuple->ndata,0,
  2268.                     disp->ntuple->ndim)] );
  2269.      xp = &(disp->ntuple->data[INDEX(0,disp->binding.x,disp->ntuple->ndim)]);
  2270.      yp = &(disp->ntuple->data[INDEX(0,disp->binding.y,disp->ntuple->ndim)]);
  2271.      nt = &(disp->ntuple->data[INDEX(0,0,disp->ntuple->ndim)]);
  2272.      
  2273.      i = 0;
  2274.      xyp = *xy;
  2275.      
  2276.      for ( ; xp<last; 
  2277.       xp += nt_dim, yp += nt_dim, nt += nt_dim )
  2278.      {
  2279.       if (disp->nt_cut!=NULL && doCuts(nt, disp->nt_cut)) continue;
  2280.       
  2281.       if (*xp >= xhigh)
  2282.       {
  2283.            if (*yp > yhigh)
  2284.             disp->bins.totals[2][2] += 1.0;
  2285.            else if (*yp < ylow)
  2286.             disp->bins.totals[2][0] += 1.0;
  2287.            else
  2288.             disp->bins.totals[2][1] += 1.0;
  2289.       }
  2290.       else if (*xp < xlow)
  2291.       {
  2292.            if (*yp >= yhigh)
  2293.             disp->bins.totals[0][2] += 1.0;
  2294.            else if (*yp < ylow)
  2295.             disp->bins.totals[0][0] += 1.0;
  2296.            else
  2297.             disp->bins.totals[0][1] += 1.0;
  2298.       }
  2299.       else
  2300.       {
  2301.            if (*yp > yhigh)
  2302.             disp->bins.totals[1][2] += 1.0;
  2303.            else if (*yp < ylow)
  2304.             disp->bins.totals[1][0] += 1.0;
  2305.            else
  2306.            {
  2307.             disp->bins.totals[1][1] += 1.0;
  2308.             
  2309.             if (!xlog)
  2310.              *xyp++ = *xp;
  2311.             else
  2312.              *xyp++ = log10(*xp);
  2313.             
  2314.             if (!ylog)
  2315.              *xyp++ = *yp;
  2316.             else
  2317.              *xyp++ = log10(*yp);
  2318.             
  2319.             if (*xp > m->xh) m->xh = *xp;
  2320.             if (*xp < m->xl) m->xl = *xp;
  2321.             
  2322.             if (*yp > m->yh) m->yh = *yp;
  2323.             if (*yp < m->yl) m->yl = *yp;
  2324.             
  2325.             i++;
  2326.            }
  2327.       }
  2328.      }
  2329.      
  2330.      return i;
  2331. }
  2332.  
  2333.  
  2334. static int drawLego2D(display disp)
  2335. {
  2336.      int rc;
  2337.      
  2338.      if (!init_done)
  2339.      {
  2340.       initPlot( disp );
  2341.       init_done = 1;
  2342.      }
  2343.           
  2344.      switch (plot_drvr)
  2345.      {
  2346.       NeXTFunc(drawLego2D,( &disp->marginRect ));
  2347.       
  2348.       PSFunc(drawLego2D,( &disp->marginRect ));
  2349.       
  2350.       UPFunc(drawLego2D,( disp ));
  2351.       
  2352.       XIVFunc(drawLego2D,(disp->bins.xAxis.nBins,
  2353.                   disp->bins.yAxis.nBins,
  2354.                   disp->bins.binMin,
  2355.                   disp->bins.binMax,
  2356.                   disp->bins.data));
  2357.       
  2358.       X11Func(drawLego2D,(disp->bins.xAxis.nBins,
  2359.                   disp->bins.yAxis.nBins,
  2360.                   disp->bins.binMin,
  2361.                   disp->bins.binMax,
  2362.                   disp->bins.data));
  2363.       
  2364.       MACFunc(drawLego2D,( &disp->marginRect ));
  2365.  
  2366.      default:
  2367.       rc = -1;
  2368.      }
  2369.      
  2370.      return rc;
  2371. }                 
  2372.  
  2373. int h_shade(display disp, float var1, float var2)
  2374. /*
  2375.  * Shade an area of a plot to indicate area. Make sure not
  2376.  * to outside the bounds of the plot.
  2377.  */
  2378. {
  2379.      int rc;
  2380.      float low, high;
  2381.      
  2382.      if (disp == NULL) return -1;
  2383.  
  2384.      low = (var1 - disp->xAxis.low)/(disp->xAxis.high - disp->xAxis.low);
  2385.      if (low < 0.0) low = 0.0;
  2386.      low = low * disp->marginRect.size.width + disp->marginRect.origin.x;
  2387.  
  2388.      high = (var2 - disp->xAxis.low)/(disp->xAxis.high - disp->xAxis.low);
  2389.      if (high > 1.0) high = 1.0;
  2390.      high = high * disp->marginRect.size.width + disp->marginRect.origin.x;
  2391.      
  2392.      switch (plot_drvr)
  2393.      {
  2394.       NeXTFunc(shade,(low,high,disp->marginRect.origin.y,
  2395.               disp->marginRect.origin.y+
  2396.               disp->marginRect.size.height));
  2397.  
  2398.       PSFunc(shade,(low,high,disp->marginRect.origin.y,
  2399.             disp->marginRect.origin.y+
  2400.             disp->marginRect.size.height));
  2401.       
  2402.       MACFunc(shade,(low,high,disp->marginRect.origin.y,
  2403.               disp->marginRect.origin.y+
  2404.               disp->marginRect.size.height));
  2405.  
  2406.       XIVFunc(shade,(low,high,disp->marginRect.origin.y,
  2407.              disp->marginRect.origin.y+
  2408.              disp->marginRect.size.height));
  2409.       
  2410.         X11Func(shade,(low,high,disp->marginRect.origin.y,
  2411.              disp->marginRect.origin.y+
  2412.              disp->marginRect.size.height));
  2413.       
  2414.      default:
  2415.       rc = -1;
  2416.      }
  2417.      
  2418.      return rc;
  2419. }                 
  2420.  
  2421.  
  2422. typedef struct fp
  2423. {
  2424.      float x,f;
  2425.      struct fp *n;
  2426. } fp;
  2427.  
  2428. static fp *pCache(int n );
  2429. typedef double (*pfun)( double, double *);
  2430.  
  2431. static int plotFunc( display disp, func_id func )
  2432. {
  2433.      int yClip = 0;
  2434.      int i, nCacheUsed;
  2435.      float ylow, yhigh;
  2436.      float maxKink;
  2437.      float delta1, minDelta;
  2438.      float delta2, kink;
  2439.      int isover;
  2440.      
  2441.      fp *list=NULL, *p=NULL, *t;
  2442.      int *over, *overp;
  2443.      float *xy, *xyp;
  2444.      
  2445. #define NSEG_START 20
  2446.  
  2447.      if (disp == NULL || func == NULL) return 0;
  2448.      if (disp->xAxis.high <= disp->xAxis.low) 
  2449.      {
  2450.       h_error("Plotfunc - xAxis.high <= xAxis.low");
  2451.       return 0;
  2452.      }
  2453.      
  2454.      /*
  2455.       * start from the beginning of the cache of point structures.
  2456.       */
  2457.      nCacheUsed = 0;
  2458.      
  2459.      maxKink = 5.0 * 3.141593/180.0;
  2460.      minDelta = 4.0 * (disp->xAxis.high - disp->xAxis.low) /
  2461.       disp->marginRect.size.width;              /* a couple of pixels */
  2462.      
  2463.      if (init_done || !disp->yAxis.flags.autoScale) yClip = 1;
  2464.      ylow = FLT_MAX;
  2465.      yhigh = -FLT_MAX;
  2466.      
  2467.      delta1 = (disp->xAxis.high - disp->xAxis.low)/((float)NSEG_START);
  2468.      for (i=0; i<=NSEG_START; i++)
  2469.      {
  2470.       if (i==0) 
  2471.       {
  2472.            list = pCache(nCacheUsed++);
  2473.            p = list;
  2474.       }
  2475.       else
  2476.       {
  2477.            p->n = pCache(nCacheUsed++);
  2478.            p = p->n;
  2479.       }
  2480.       p->x = disp->xAxis.low + ((float) i) * delta1;
  2481.       p->f = (*(pfun)(func->funcPtr))((double)p->x,
  2482.                       (double *)func->paramBlk);
  2483.       p->n = NULL;
  2484.      }
  2485.      
  2486.      for (p=list; (p->n)->n != NULL; p = p->n )
  2487.      {
  2488.       kink = atan((((p->n)->n)->f - (p->n)->f)/
  2489.               (((p->n)->n)->x - (p->n)->x)) -
  2490.                atan(((p->n)->f - p->f)/((p->n)->x - p->x));
  2491.       delta2 = ((p->n)->n)->x - (p->n)->x;
  2492.       delta1 = (p->n)->x - p->x;
  2493.       while ( fabs(kink) > maxKink && 
  2494.          (delta1 > minDelta || delta2 > minDelta) )
  2495.       {
  2496.            if (delta1 <= delta2)
  2497.            {
  2498.             t = pCache(nCacheUsed++);
  2499.             t->x = (p->n)->x + delta2/2.0;
  2500.             t->f = (float) 
  2501.              (*(pfun)(func->funcPtr))((double)t->x,
  2502.                         (double *)func->paramBlk);
  2503.             t->n = (p->n)->n;
  2504.             (p->n)->n = t;
  2505.            }            
  2506.            if (delta1 >= delta2)
  2507.            {
  2508.             t = pCache(nCacheUsed++);
  2509.             t->x = p->x + delta1/2.0;
  2510.             t->f = (float) 
  2511.              (*(pfun)(func->funcPtr))((double)t->x,
  2512.                         (double *)func->paramBlk);
  2513.             t->n = p->n;
  2514.             p->n = t;
  2515.            }            
  2516.            kink = atan((((p->n)->n)->f - (p->n)->f)/
  2517.                (((p->n)->n)->x - (p->n)->x)) -
  2518.                 atan(((p->n)->f - p->f)/((p->n)->x - p->x));
  2519.            delta2 = ((p->n)->n)->x - (p->n)->x;
  2520.            delta1 = (p->n)->x - p->x;
  2521.       }    
  2522.      }
  2523.      
  2524.      /* 
  2525.       * We now have the complete list of points to plot. Create the
  2526.       *  list of x-y points and overflow flags.
  2527.       */
  2528.      xy = (float *) malloc( 2 * nCacheUsed * sizeof(float) );
  2529.      over = (int *) malloc( (nCacheUsed+1) * sizeof(int) );
  2530.      xyp = xy;
  2531.      overp = over;
  2532.      
  2533.      for (p=list, i=0; i<nCacheUsed; p = p->n, i++ )
  2534.      {
  2535.       isover = 0;
  2536.       
  2537.       if (disp->xAxis.flags.log)
  2538.            if (p->x > 0.0)
  2539.             *xyp++ = p->x;
  2540.            else
  2541.            {
  2542.             isover = 1;
  2543.             *xyp++ = -FLT_MAX;
  2544.            }
  2545.       else
  2546.            *xyp++ = p->x;
  2547.  
  2548.       if (disp->yAxis.flags.log)
  2549.            if (p->f > 0.0)
  2550.             *xyp++ = p->f;
  2551.            else
  2552.            {
  2553.             isover = 1;
  2554.             *xyp++ = -FLT_MAX;
  2555.            }
  2556.       else
  2557.            *xyp++ = p->f;
  2558.       
  2559.       if (isover || 
  2560.           (yClip && (p->f > disp->yAxis.high || p->f < disp->yAxis.low)))
  2561.            *overp++ = i;
  2562.       else
  2563.       {
  2564.            if (p->f > yhigh) yhigh = p->f;
  2565.            if (p->f < ylow && (!disp->yAxis.flags.log || p->f > 0.0))
  2566.             ylow = p->f;
  2567.       }
  2568.      }
  2569.      *overp++ = i;
  2570.      
  2571.      /*
  2572.       * Set scales if it has not been done already.
  2573.       */
  2574.      if (! init_done)
  2575.      {
  2576.       if (disp->yAxis.flags.autoScale)
  2577.       {
  2578.            disp->yAxis.low = ylow;
  2579.            disp->yAxis.high = yhigh;
  2580.            h_adjustAxis(&(disp->yAxis.low), &(disp->yAxis.high),
  2581.                 0, disp->yAxis.flags.log );
  2582.       }
  2583.       initPlot( disp );
  2584.       init_done = 1;
  2585.      }
  2586.      
  2587.      /* 
  2588.       * plot the line, clean up and we're done!
  2589.       */
  2590.      drawLine1D( disp, xy, over, nCacheUsed, 0, func->lineStyle );
  2591.      
  2592.      free(xy);
  2593.      free(over);
  2594.      
  2595.      return 0;
  2596. }
  2597.  
  2598. static fp *pCache(int n )
  2599. {
  2600.      static npoints = 0;
  2601. #define FP_BLK 100
  2602.      static fp *list;
  2603.      
  2604.      if (npoints == 0)
  2605.      {
  2606.       if ( (list = (fp *)malloc( FP_BLK * sizeof(fp) )) == NULL)
  2607.       {
  2608.            h_error("Plotfunc is out of memory");
  2609.            return NULL;
  2610.       }
  2611.       npoints = FP_BLK;
  2612.      }
  2613.      
  2614.      if (n >= npoints)
  2615.      {
  2616.       if ( (list = (fp *)realloc( list, 
  2617.                      (npoints+FP_BLK) * sizeof(fp) )) == NULL)
  2618.       {
  2619.            h_error("Plotfunc is out of memory");
  2620.            return NULL;
  2621.       }
  2622.       npoints += FP_BLK;
  2623.      }
  2624.      
  2625.      return &list[n];
  2626. }
  2627.